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We extend the imaginary-time formulation of the equilibrium quantum many-body theory to 
steady-state nonequilibrium with an application to strongly correlated transport. By introducing 
Matsubara voltage, we keep the finite chemical potential shifts in the Fermi-Dirac function, in 
agreement with the Keldysh formulation. The formulation is applied to strongly correlated transport 
in the Kondo regime using the quantum Monte Carlo method. 
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A coherent formulation of equilibrium and nonequilib- 
rium is one of the ultimate goals of statistical physics. In 
the last two decades, this has become a particularly press- 
ing issue with the advances in nanoelectronics. Although 
it has long been considered such Gibbsian description 
may exist in the steady-state nonequilibrium imple- 
mentation of time-independent nonequilibrium quantum 
statistics has produced limited success [2] without widely 
applicable algorithms. 

In nanoelectronics, the strong interplay between many- 
body interactions and nonequilibrium demands nonper- 
turbative treatments of the quantum many-body effects. 
Perturbative Green function techniques ^ 3] have been 
successful, but are often plagued by complicated dia- 
grammatic rules and are limited to simple models. In 
the last few years, important advances have been made in 
this field to complement the diagrammatic theory. Time- 
dependent renormalization group ^ and density- 
matrix renormalization group method [7| were applied to 
calculate the real-time convergence toward the steady- 
state. Real-time methods 0, 01 calculate the process 
toward the steady-state and therefore have clear physi- 
cal interpretations. Unfortunately they often suffer from 
long-time behaviors associated with low energy strongly 
correlated states and finite size effects. Direct construc- 
tion of nonequilibrium ensembles through the scatter- 
ing state formalism 0, [^, [l^ and field theoretic ap- 
proach [11] have provided new perspectives to the prob- 
lem. 

The main goal of this work is to provide a critical 
step toward the time-independent description of equi- 
librium and steady-state nonequilibrium quantum statis- 
tics. In addition to the resolution of this fundamental 
problem, we provide a strong application. The steady- 
state nonequilibrium can be solved within the same for- 
mal structure as equilibrium, and therefore the power- 
ful equilibrium many-body tools, such as the quantum 
Monte Carlo (QMC) method, can be easily applied to 
complex transport systems with many competing inter- 
actions. We demonstrate this point by applying this 
formalism to strongly correlated transport in the Kondo 
regime by using QMC. In contrast to the real-time meth- 



ods, this approach starts from the steady-state and sim- 
ulates the effect of many-body interaction. However, nu- 
merical analytic continuation and low temperature cal- 
culation, especially with the QMC application, are tech- 
nical difficulties. 

In the following, we first construct a time-independent 
statistical ensemble of steady-state nonequilibrium [2] in 
the non-interacting limit with the introduction of Mat- 
subara voltage. We show that the interacting imaginary- 
time Green function can be mapped to the retarded 
Green function after an analytic continuation to the real- 
time and real-bias. The spectral representation is used 
to carry out the numerical analytic continuation. We 
use QMC for the Kondo dot system [13] to solve for the 
strongly correlated transport. 

The expectation value of an operator A is defined on 
the ensemble propagated from the remote past. 



(A) = Hm 



Tr[p(T)i] 
Trp(r) ■ 



(1) 



with p{T) — e'^^'^ poe~'^^^ where the initial non- 
interacting ensemble in the remote past is given by po- 
The total Hamiltonian is given hy H = Ho + V with the 
non-interacting part 



Hn = 



h.C.) 



(2) 



where cj^^,^ is the conduction electron creation operator 
on the a reservoir (a = 1 for the source and a = — 1 for 
the drain leads) with the continuum index k and spin a. 

It is crucial that we choose the initial ensemble to be 
a fully established steady-state nonequilibrium. Since we 
consider an open system with infinite volume, the time- 
evolution of a zero-current ensemble after any finite time 
t, however long, retains the non-vanishing contribution 
from the remote past, as pointed out by Duyon and An- 
drei i. 

For the moment, let us consider the noninteracting 
model Hq. The time-evolution of the nonequilibirium 



2 



steady-state ensemble is given by Hershfield ^, ^ with 

po = e-/5(^o-«>o)^ (3) 

where the operator Yq imposes the nonequihbrium 
boundary condition in terms of the scattering states of 
Hq. In the non- interacting system the scattering states 



V-'ifeo- calculated explicitf 

Lippmann-Schwinger equation [13 




in the form of the 



a' k'(T 



(4) 



where gdi^) is the retarded Green function of the quan- 
tum dot (QD) site. For an infinite band system, 
(/ci(e) becomes = (e — + with the hy- 

bridization broadening F = Fl + Tfj, where Fq, = 
7ri^7V(0) [Af(0)=density of states of the leads]. It can 
be shown in a straightforward calculation that Hq = 
J2aka^aki>ika'^aka- The boundary condition operator 
Yq imposes the nonequilibrium by shifting the chemical 
potentials to the scattering states V'lfco- i^'^^ ^^'^ bare 



conduction electrons cl 



with 



Yn 



(5) 



We have chosen the voltage drop to be symmetric about 
the QD region, although by choosing 7^ we can apply 
the following formalism in general. 

The expectation value {A), Eq. ([1]), is expressed as 

(A) ^ (^J I?[V't,V']^(V'UO),V'(0))e'-''^^*)''*) , where the 
average is performed with respect to po- The Lagrangian 
is L{t) = Eafe^V'L^(*)(«^t - ^akHakait). By defining 



q:$/2, we have po 



andL(i) = J2c.ka ^lkAt)i^dt-iak-a<^/2yj^kAt)- Note 
that the states on the Fermi energy in each lead {cak — 0) 
have the different time-evolution rates, 

In order for the analytic continuation to work, the ex- 
tra time-evolution rate is factored out formally as 



Ipakcrit) = e 



'ipakait), 



(6) 



which does not affect po, but changes the Lagrangian to 

= J2aka '^lkJt)i^^t - eafc)V'afc^(i)- 

Now we introduce the analytic continuation with it 
r for the field variables 'ipakait) and ipl,j^„{t). The crucial 
step is to realize that the phase factor in Eq. ^ becomes 
divergent (or vanishing) in e~"*'^/^ and that this can be 
avoided by introducing the Matsubara voltage^ 



. 47rm 
iip^ ^ $ with ip^ = —— (m 

P 



integer). (7) 



The bosonic Matsubara frequency guarantees the same 
periodic boundary condition of thermal Green functions 



as the equilibrium formalism. Here we have two analytic 
continuations, one in time and the other in bias. Fend- 
ley et al [31 has first introduced the Matsubara voltage 
for the bare reservoir states within the Bethe Ansatz for- 
malism. However, when implemented in Green function 
theory [l^ discrepancies from the Keldysh method have 
been pointed out. 

The time-ordered QD Green function is defined as 
Sddi''') = ^('^'^(''')'^^(O)) where the propagation in 
the imaginary-time is given by the action So{t) — 

'Eaka^ikA'^)i9r " Cafe ) V'afe^r (t) = Eaferr V'Lrr M I^r - 

^ak — %{i^m — ^)]'^ak(j{T)- Here, the evolution in 
the imaginary-time is governed by the effective non- 
interacting Hamiltonian Kq = Hq + {iiprn. — ^)yo- Using 
the expansion of the scattering states [10|, the Fourier 
transformation of S2(;(*w„) at the Matsubara frequency 
Un = (2n 



l)7r//3 can be readily calculated as 
Fa/F 



(8) 



with Tnm — F • sign(Li;„ — aLpm/^). With the analytic 
continuations i^Pm — *■ ^ followed by iojn lu + irj, we 
recover the retarded Green function gd(j^)- 

With an interaction V, the effective action is S* = 
S'o — drV [dj. (r) , da- (r)] or equivalently the effective 
Hamiltonian K becomes 



K = Ko + V = Hq + {i(p„, - $)ro + V. 



(9) 



Now we show that the imaginary-time evolution through 
K followed by the analytic continuation itpm and 
then iLUn — > LU+irj gives the same retarded Green function 
calculated in real time. 

The imaginary-time Green function is expanded in the 
interaction picture with 6"^^° as the unperturbed den- 
sity matrix and V/(r) = e''^»Ve^^°, 



Qddir) 



Tr 


PQ%e-fo^ 


^'^^(^')d(r)rft(0) 


Tr 


PQ%-e 


^ f^^ dr'Vr(T') 





(10) 



First consider the first order expansion of the numera- 
tor for T > 0. We denote the imaginary-time ordering 
of a followed by b as (6a)/. Thus there are two possible 
first order contributions, (dVd^)i and {Vdd'')]. An ex- 
plicit expression for {dVd^)i after Fourier transformation 
becomes 

idVd^)i = X(nM|m) J"^'^'^^ {k\d^\n)x (11) 



KQm — Kok 
POn 



POn 



ILOn + KQn - Kq^ iiOn + Kqu - Kqu 
POm POk 



iuJn + Kon - Kom iuJn + Kon - Ki 



Ok 
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with respect to the unperturbed energy eigenstate \n) at 
the eigenvalue K^n = E^n + {i'Pm — ^)Yon- In the above 
derivation we used the critical relation 



-/3[ifo + (iVm-*)yo] 



Po, 



(12) 



which holds only when (pm is a Matsubara frequency 
and [Ho,Yo] = 0. Since \n) can be constructed from 
the scattering states, the eigenvalues Ybn are half-integer 
with e-''''^™^"" = 1, hence e-''I-^0"+(*^--*)^«"l = 

Now we consider the real-time retarded Green function 
GUt)^m[G>{t)-G<it)] with 



in the algorithm are the initial Green function Eq. ^ 
and multiple runs performed at different ipm- In the 
QMC calculations, the discretization error (FAt = 0.2) 
makes high frequency quantities unreliable and we thus 
have limited ipm up to 1.5U. Throughout this paper, 
the unit of energy is given by the hybridization strength 
T = Tl+Tr = 1. 

We start the numerical analytic continuation by study- 
ing the analytic structure of the self-energy in the second 
order at {iujn,iipm) 



G>{t) 



Tr 



PoTKe-*^K'i*'^Ht')d(i_)rft(o+) 



n 

.1=1 



(15) 



TV poTKe-'fKdt'Vj(t') 



(13) 



/c[i(l f 02) fas + (1 /qi)/q2(I faa) 



{ai - 0:2 + as) 



ei + 62 - €3 



where the time-ordering is defined on the Keldysh con- 
tour K. 0+ is on the first half of the contour K, (— T 
T), and t- is on (T -T). G<{t) is similarly defined. 
In the first order, we have 6 distinct time-ordering in 
G^{t) along K, namely (ddty)^, idVd^)K, iVdd^)K, 
{d^dV)K, {d^Vd)K and {Vd^d)K- An exphcit calcula- 
tion for {dVd^)K after Fourier transformation gives 

{dVd^)K = V(nM|m)-M5^(fc|dt|n) X 



nrnk 



(14 



POri 



POn 



Lo + Eon - Eo,n + irj LU + Eon ~ Eok + iri 

This expression agrees with the first two terms in the 
parenthesis in Eq. pip after the analytic continuations 
ifm - $ ^ (i.e. - i^om - £^om etc.) and 

iojn ^ oj + irj. Similarly, {Vd'^d)K^ a cyclic permutation 
of {dVd'^)i, produces the third term in Eq. pT|) and an- 
other cyclic permutation {d'^dV)K gives the last term. 
The remaining real-time-orderings {d'^dV)Ki {d'^Vd)K 
and {Vd'^d)K are generated by the cyclic permutations 
of the imaginary-time ordering {Vdd'^)i. Such a mapping 
can be established in the higher order expansions. For 
instance, in the second order of V , the 3 distinct order- 
ings {VdVd'^)!, (yVdd^)i and {dVVd'<)i produce the 12 
distinct real-time orderings {VdVd^)K, {Vd^Vd)K, etc. 

The above mapping between the real- and imaginary- 
time Green functions is expected since the term-by-term 
correspondence remains the same regardless of the val- 
ues of iifm — ^ and the equilibrium limit guarantees 
the equivalence of perturbation expansion in both ap- 
proaches. The main effect of the Hamiltonian Eq. ([9]) is to 
correctly give the initial statistics by Hq ~ ^Yq [Eq. 
and the time-evolution by H after iipm — * [17]. 

From now on, we discuss the numerical implementa- 
tion of the above formulation to the Kondo anomaly us- 
ing the QMC method. In this work, the Hirsch-Fye [l3| 
algorithm is applied to the on-site Coulomb interaction 
Hi = U {ud] — 5) {ndi — W The only modifications 



Here /^^ = /(e^ — aif-), the Fermi-Dirac function with 
the shifted chemical potential. This expression can be 
derived with the standard equilibrium second order per- 
turbation theory [l9| but with the nonequilibrium Green 
function Eq. ([8]) as an input. Similarly to Eq. (fT2|) . the 
critical step /(e + a "^"^"'^ ) = /(e — af-) has been used. 
After taking i^Pm ~* ^ and then iujn lo + irj, this ex- 
pression maps to the correct retarded self-energy in the 
Keldysh formalism [20]. 

' Motivated by the form of the above self-energy, we 
decompose the numerical self-energy in a spectral repre- 
sentation with multiple branch-cuts with respect to e. 



de- 



cr^(e) 



7- 



(16) 



with odd integers 7. We fit the spectral function (J^{e) 
defined on a logarithmic frequency mesh. In the fit we 
used I7I < 9 {i.e. 8 branch-cuts ^). FIG. [Hb) shows 
the analytic continuation (iipm ^) of the perturbation 
self-energy Eq- after the fit has been found. 

After benchmarking the analytic continuation, we an- 
alyze the QMC self-energies shown in FIG. [T]^c) for ipm 
with m — 0, - ■ ■ ,7. In FlCf^d), the spectral function for 
G-^'{u!) is plotted with bias <& from to 0.5 with an inter- 
val 0.05. The equilibrium Kondo peak becomes quickly 
quenched at $ ^ Tk with Tk ~ 0.11 determined at 
HWHM for $ = [ttTA^Tk) = 0.5]. As $ increases 
further, side peaks develop in agreement with the fourth 
order perturbation results [13] . The position of the peaks 
roughly scales with $ and their width is comparable to 
r, indicating that they have weak correlation effects. 

With the spectral function for G^{oj), one can calcu- 
late the current from the relation [2l| 

^ - ^ / de[G^ie) - G^ie)] [/^(e) - fn{e)] . (17) 

FIG. shows the differential conductance as a function 
of <I>. The thin solid line is the non- interacting limit. 
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FIG. 1: (a) Second order perturbation self-energy at 
{iuj„,iifrn) with Matsubara voltage ipm = 4m7r//3. (b) An- 
alytic continuation iipm ^ ^ performed based on the fit, 
Eq. (|16|) . The result agrees very well with the exact continu- 
ation from the analytic expression, (c) Self-energy calculated 
from the quantum Monte Carlo method, (d) Spectral func- 
tion A{u!) of the QD Green function. The Kondo peak quickly 
disappears as the bias "1> is increased and develops into two 
broad peaks at cij ~ ±$. 



The $ = U/2 peak corresponds to the co-tunnehng with 
the charge-excited QD. The $ = /7 peak is due to the 
inelastic QD-lead tunneling. 

We have plotted the zero-bias ($ = 0, ?7 = 5) spec- 
tral function (dashed line), -Aeq(u;), calculated from the 
maximum entropy method [23| . As expected, the finite- 
bias peak width is much narrower than the equilibrium 
prediction due to the destruction of the Kondo peak at 
finite dc-bias. For comparison, the frequency lu is scaled 
to $/2 to match the chemical potential profile of the 
source-drain with respect to the QD. 

With increased U = 10, the anomalous Kondo peak 
becomes sharper with the HWHM ^hwhm <C Tk for 
Tk estimated from the spectral function in FIG. [2l In 
addition to the Kondo and inelastic charge peaks, the 
transport across the side peaks in FIG. [Hd) emerges as 
a weak peak between $ = and $ = U/2. 

We acknowledge support from the National Science 
Foundation DMR-0426826 and computing resources at 
CCR of SUNY Buffalo. 
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FIG. 2: DC-conductance of Kondo quantum dot system. The 
width of the anomalous peak is significantly narrower than 
what the zero-bias spectral function predicts [Acqico,^ = 0) 
with u scaled according to $/2, short-dashed line], due to 
the destruction of Kondo resonance at finite bias In the 
strongly correlated regime {U — 10), the Kondo peak becomes 
more pronounced with strong temperature-dependence. At 
$ ~ f//2 and U, the broad inelastic transport peak emerges. 
The non-interacting limit {U = 0) is shown as thin line. 



[/ = 0. With the chemical potentials displaced by ±<I>/2 
from the QD level, the HWHM occurs at $/2 w F = 1. 

As the interaction is turned on, the zero-bias conduc- 
tance becomes narrower. At J7 = 5 (solid circle), the 
anomalous Kondo peak begins to develop. The zero-bias 
limit approaches the unitary limit as T — > 0. At higher $, 
inelastic transport peaks appear at $ = J7/2 and $ = [/. 
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